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The radiative recombination rates of interacting electron-hole pairs in a quantum dot are strongly 
affected by quantum correlations among electrons and holes in the dot. Recent measurements of the 
biexciton recombination rate in single self-assembled quantum dots have found values spanning from 
two times the single exciton recombination rate to values well below the exciton decay rate. In this 
paper, a Feynman path-integral formulation is developed to calculate recombination rates including 
thermal and many-body effects. Using real-space Monte Carlo integration, the path-integral expres- 
sions for realistic three-dimensional models of InGaAs/GaAs, CdSe/ZnSe, and InP/InGaP dots are 
evaluated, including anisotropic effective masses. Depending on size, radiative rates of typical dots 
lie in the regime between strong and intermediate confinement. The results compare favorably to 
recent experiments and calculations on related dot systems. Configuration interaction calculations 
using uncorrelated basis sets are found to be severely limited in calculating decay rates. 

PACS numbers: 78.67.Hc,31.15.Kb 



I. INTRODUCTION 

The small size and strong optical properties of self- 
assembled quantum dots (QDs) make them appealing 
candidates for optoelectronic devicesi^ When light is 
absorbed, photons create electron-hole (eh) pairs (exci- 
tons) that become confined in the quantum dot. Re- 
cent photoluminescence (PL) spectra have measured the 
recombination energy of electron-hole pairs with mcV 
resolution!^ Analysis of single dot PL spectra at dif- 
ferent incident light intensities reveals that the exciton 
recombination energy is shifted by other "spectator" ex- 
citons and free charges in the dotA^ For example the 
recombination energy is red-shifted a few meV by the 
presence of a spectator excitoni^ Detailed understand- 
ing of the effect of spectators on recombination is impor- 
tant for non-linear optical applications, such as quantum 
logic gates^ or turnstilesAl 

The rates of the PL processes determine the steady- 
state occupation of the dots for a given incident 
intensity^ Time-resolved photoluminescence measure- 
ments can track the electron-hole recombination rate in 
single self-assembled quantum dots^ Recent experiments 
give differing results about the decay rate of the biexci- 
ton relative to that of an isolated exciton in the same dot. 
Measurements on a single CdSe/ZnSe dot find a biexci- 
ton decay rate Txx about equal to the exciton rate Tx r 
while other experiments on similar sized CdSe/ZnSe QDs 
report a biexciton decay rate twice the exciton rate£ Sim- 
ilar measurements in InGaAs have found Txx /Tx ~ 1 .5^ 
Txx Ax -2,12*11 and even T xx /T x ^0.33i2, 

Theoretically, there are two limits to consider for re- 
combination rates. In the strong confinement limit, the 
exciton and biexciton wave function is a simple product of 
the electron and hole single-particle wave functions in the 
dot. Coulomb interactions are assumed to only slightly 
perturb the wave function. In that case the recombina- 



tion rates contain matrix elements of the single-particle 
wave functions, which are the same for excitons and biex- 
citons. Taking into account the number of allowed decay 
channels, the biexciton should decay at twice the exci- 
ton rate, Txx/Tx = 2. The other limit is the weak 
confinement limit, which applies when the exciton bind- 
ing energy significantly exceeds the single-particle level 
spacing of the dot. In that limit, the exciton or biexciton 
is bulk-like, bound together as a small composite parti- 
cle. This exciton or biexciton unit is weakly confined in 
a dot much larger than the exciton or biexciton radius. 
In this case the dipole matrix element is dominated by 
the exciton or biexciton structure, which is independent 
of dot size. The composite particle has a coherent wave 
function that extends across the volume of the dot, lead- 
ing to constructive addition of radiative matrix elements 
for exciton decay. Thus, in the weak confinement limit, 
the radiative decay rate of the exciton increases with dot 
size, until the dot diameter approaches the wavelength 
of the emitted light. For the biexciton, the exciton final 
state after recombination suppresses this constructive en- 
hancement, significantly reducing the value of Txx /Tx 
in the weak confining limit. In the intermediate regime, 
the exciton wave function generally cannot be separated, 
except for some special choice of the external potential, 
such as a harmonic confinement^ Still, the coherent ex- 
tent of the many-particle wave function — the coherence 
volume*^ — leads to an increasing decay rate with increas- 
ing dot size and will play an important role in the inter- 
pretation of our results. In this paper we show that the 
radiative decay rates of typical self-assembled dots lie in 
the regime between strong and intermediate confinement. 

Theoretical descriptions of single particle (electron or 
hole) states in quantum dots have improved greatly in 
the last ten years ^i^i^i^ y 6 t the description of exci- 
ton or multi-exciton states is not as well developed. The 
energies of states with several electrons and holes are 
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usually treated within first-order perturbation theory. 
Some spectral energies, such as the biexciton shift, re- 
quire treatment of correlation with configuration interac- 
tion (CI) or quantum Monte Carlo (QMC) techniques^ 
The limited accuracy of approximate CI wave functions 
is known to affect the calculated energies J£ 

There have been a few attempts to calculate biexci- 
ton decay rates in quantum dots. Takagahara used a 
variational calculation to determine the decay rate of 
a biexciton in an infinite barrier spherical dot with di- 
electric effects^ More recently, Ungier et al. have calcu- 
lated the biexciton decay rate for zincblende and wurtzite 
structures iiS Using CI expansions, Corni et al. have stud- 
ied the size dependence of exciton and biexciton recom- 
bination rates in strain-induced dotsiSi These dots are 
formed in a near-surface InGaAs/GaAs quantum well by 
the stress field of an InP self-assembled island grown on 
the surface. These dots have much shallower confinement 
than self-assembled dots, are often much larger, and are 
well-approximated by truncated 2-d parabolic confine- 
ment. This puts the strain induced dots well into the 
weak confinement regime, contrary to the more common 
self- assembled dots that are subject of this paper. Re- 
cently, Narvaez et al. have performed CI calculations on 
InGaAs/GaAs self-assembled dots beyond the effective 
mass approximation using pscudopotentials^i These CI 
results must be viewed with some caution since decay 
rates are more sensitive than energies to errors in the 
wave function, as we will show in this paper. 

In this paper we develop a Fcynman path integral de- 
scription of exciton and biexciton recombination rates. 
This technique can be easily applied to complicated dot 
geometries, does not depend on a finite basis set, and 
fully treats correlation. In Sec.|n]we derive a path inte- 
gral expression for the recombination rate. This expres- 
sion is then evaluated using a real-space Monte Carlo 
technique that we introduce in Sec. IIIII In Sec. IIVI we 
apply our path integral technique to a model system and 
compare with full CI calculations. In Sec. |V] we apply 
both the path integral technique and the CI expansion 
to realistic three-dimensional models of InGaAs/GaAs, 
CdSc/ZnSc, and InP/InGaP dots and compare to sin- 
gle dot experiments. While our path integral method is 
currently restricted to single-band effective mass approx- 
imation (EMA) models, the insights, trends, and even 
quantitative rates revealed in these make them quite use- 
ful, as we conclude in Sec. IVII 



II. METHOD 

Starting from a standard treatment of electron-hole 
radiative recombination in the effective mass approxima- 
tion, we rewrite the square of the matrix element in the 
rate equation as a path integral expression. In the path 
integral formalism, we will show that the rate is propor- 
tional to the ratio of two path integrals: one with the 
standard thermal trace, and the other with a "radiating" 
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FIG. 1: (Color online) Illustration of our path-integral cal- 
culation of biexciton recombination rates in self-assembled 
quantum dots. We express the rate as a ratio of path in- 
tegrals with (a) diagonal and (b) radiating constraints, see 
Eq. ©. We evaluate the path integrals using Monte Carlo 
integration on realistic three-dimensional models. A typical 
path contributing to the integrals for an InGaAs/GaAs dot 
is shown in (c). These paths sample the probability density, 
energy, rate, and other properties of the radiating states. 



configuration that pairs an electron and hole, as illus- 
trated in Fig. Ofa) and (b). 



A. Exciton recombination rate within the effective 
mass approximation 

Our Hamiltonian is a commonly used effective mass 
model, 

where V e and Vh describe a lens-shaped confining po- 
tential and the hole effective masses mn are anisotropic. 
In contrast to previous approaches to this problem, we 
do not construct single-particle or variational wave func- 
tions. Rather, we use Metropolis Monte Carlo to sam- 
ple the recombination rate directly from a path integral. 
The path integral Monte Carlo (PIMC) method allows 
us to calculate the density matrix for the Hamiltonian, 
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Eq. JIJ. As we describe below, this is an essentially ex- 
act solution without basis set problems or the difficulties 
of variational approaches. A snapshot of a typical path 
for an electron-hole pair in our simulation is shown in 
Fig. ^c). A sum over all such paths is a complete quan- 
tum mechanical solution for the model, Eq. |T]l. 

The rate of spontaneous decay of an exciton into a pho- 
ton is the sum of the rates of all possible decay processes. 
For generality, wc consider a state 3>f with N electron- 
hole pairs decaying to a state <I>y with N — 1 pairs. The 
rate of spontaneous emission into a photon with polar- 
ization A, momentum Kk, and energy So;, in a medium 
with index of refraction n (« \/e), is 
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which follows from Fermi's golden rule. Since the emitted 
photon has energy slightly less than the band gap (typ- 
ically 1-3 eV), the photon wavelength is much greater 
than the dimensions of the self-assembled dot (typically 
5-50 nm), so we take the k — *• limit in the current op- 
erator, jj^g y The usual approximation for the exciton 
decay rates in semiconductors is to use the envelope ap- 
proximation, in which the single particle wave functions 
are approximated as an envelope times a periodic Bloch 
function, <p(r) = ip(r)u(r). Then the current operator 
splits into a delta function on the envelope and the cur- 
rent operator j = p/m on the Bloch function. The mo- 
mentum matrix element between the conduction band 
(CB) and valence band (VB) Bloch functions is given 
by the Kane parameter, Ep = | (CB|p|VB) \ 2 /2m. Since 
all significant transitions occur in an energy range given 
by the coulomb interaction (a few tens of meV), which 
is much smaller than the gap energy, we take the usual 
approximation hw E gap . Thus, within the envelope 
approximation, the recombination rate due to transition 
a is approximately 



2nE gap Epe 
3h 2 c 3 m 
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where the point contact matrix element I is the overlap 
integral of the initial and final envelope functions, 
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B. Path integral expression for the rate 

The determination of exciton and biexciton recombi- 
nation rates using Eqs. J3J) and (0} faces two difficulties. 
First, the initial and final states contain several inter- 
acting particles, for which correlation must be treated 
carefully. Second, the total rate is a sum over all pos- 
sible transitions, T = ^2 a T a . These difficulties may be 
treated explicitly for model systems (such as harmonic 



oscillators) , but generally the direct determination of the 
matrix elements and rates in a wave function representa- 
tion leads to approximations. For example, Takagahara's 
variational calculation^ treats correlation very well for a 
single transition from a biexciton to the exciton ground 
state, but is limited to very symmetric, spherical QDs. 
On the other hand, Ungicr et al*£> treat structural de- 
tails of the dot with much care (even beyond the en- 
velope function above), but the correlation is only par- 
tially included, an approximation known to underesti- 
mate biexciton binding energies^ CI calculations can in 
principle solve the full many-particle Schrodinger equa- 
tion for excitons and biexcitons for EMA models^ and 
pscudopotcntials^I but may be severely limited by the 
underlying basis set, as we will show. 

We now derive a path-integral solution for the total 
recombination rate that will fully treat correlation, ther- 
mally distribute the initial states, and include all final 
states. To relate Eqs. © and to a path integral, we 
begin by squaring the point contact matrix element, 



l^l 2 = 



J J P%(R-N, R-Jv)PAT-l(R-iV-l, R-iV-l) 



5(r e N - r N )6(v e N - v h N )dR N dli'j 



(5) 
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where p N and p N _i are the density matrices of the ini- 
tial and final states. As in Ref. 0, we assume that 
the carriers reach thermal equilibrium before the tran- 
sition, and use the thermal density matrix of N electron- 
hole pairs, pn(Rn ,R' N ; (3). The final state can take 
on any value, so we sum over all final states, yield- 
in g/9J v-i(RAr_i,RiV-i) =^ JV - 1 )(Riv-i-R A r_ 1 ). Af- 
ter integrating out the R^y coordinates in Eq. JSJ) with 
this delta function and using Eq. we find the 

temperature-dependent radiative recombination rate, 



Tn(J3) 
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where 



(\In\ 2 )p — ^^JJ Pn(Rn, R'at; 0)S(R' n _ 1 — Rjv_i) 
S(r N - r N )5(r c N - r h N )dR N dR' N . 

(7) 

In this equation Z^ = Trpjy is the partition function for 
N electron- hole pairs and is needed to normalize pn in 
the integral. 

The thermal density matrix in Eq. 10 may be repre- 
sented as a real-space Feynman path integral^ 



p(Rn, R'at; 0) = JvR N (t)cxp -^J^Hdt 



(8) 



where the ends of the paths are Rtv(0) = R'^ and 
Rjv03) = Rat- Thus the partition function Z^ and 
the recombination integral (\In\ 2 )p can be represented 
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by path integrals that differ only by constraints on the 
paths, 



Z N = J VR N (t)cxp 

diagonal 

Z N (\I N \ 2 )p = J VR N (t)cxp 

radiating 
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(9) 
(10) 



III. COMPUTATIONAL METHODOLOGY 

The path integral expression for the recombination rate 
can be directly sampled with Monte Carlo integration, for 
two- or four-particle interacting quantum systems. We 
have implemented this as a computer simulation that al- 
lows for anisotropic masses and any three-dimensional 
confining potential we choose. 



A. Path integral Monte Carlo 



The diagonal constraint is the usual trace, Rjv(O) = 
Rjv(/3), illustrated in Fig. da). The radiating con- 
straint is a trace over the non-radiating pairs, Rjv-i (0) = 
Rat-i(/?), and a pairing of the recombining particles, 
r c N = and = r^f, as illustrated in Fig.^b). 
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It is insightful to consider how this path integral for- 
malism for the recombination rate, Eq. ©, relates to 
the strong and weak confinement limits. Consider the 
t = slice in imaginary time. In the diagonal bound- 
ary conditions, the path integral samples the diagonal 
of the density matrix in the position basis. For a non- 
interacting exciton or biexciton, the electron and hole 
sample the probability density functions of the single par- 
ticle electron and hole ground states. In the radiating 
boundary condition, the electron and hole are forced to 
coincide, but may sample two different points for t = 0_ 
and t = 0+. The effects approximately cancel out, giv- 
ing (|/jv| 2 )/3 ~ 1) appropriate for the strong confinement 
limit. With the attractive eh-interaction in the weak con- 
finement limit, the electron and hole pair together in an 
exciton. In the diagonal boundary conditions, the volume 
sampled by the electron and hole is the dot volume Vdot 
times the exciton volume ~ a x . In the radiating bound- 
ary conditions, the volume sampled is Vdot for t = 0_ 
times another factor of Vdot for t = 0+. This gives 
(\In\ 2 )[3 ~ Vdot /ax, appropriate for the weak confine- 
ment limit with dot diameter much less than the wave- 
length of light. 

Now consider a bound biexciton. One exciton has ra- 
diating boundary conditions and the other exciton has 
diagonal boundary conditions. For the strong confining 
case we see a similar cancellation of boundary condition 
effects as for the single exciton. Since we have contribu- 
tions from pairing either the spin-up or spin-down elec- 
trons and holes, we see r^x /Yx ~ 2. In the weak confin- 
ing limit, the ch-pair in the radiating boundary condition 
is bound to the other eh-pair in the diagonal boundary 
condition, with a biexciton radius axx- This binding 
suppresses a factor of Vdot in the biexciton rate, leading 
to a reduced relative rate, Yxx /Yx ~ 2a^- x /Vdot- While 
this ratio may drop below one for very large dots, most 
self- assembled dots are not much bigger than biexcitons, 
so we would not expect to see this limit except in extreme 
cases. 



With the use of Monte Carlo integration, the path inte- 
gral approach allows an essentially exact numerical solu- 
tion to many quantum statistical problems^ Quantum 
Monte Carlo methods have been useful for problems re- 
lated to this one, such as trion binding energies in quan- 
tum wellsjSiS multi-exciton energies in quantum dotsii 
and positron-electron annihilation rates^ as well as bulk 
phenomena, such as exciton-exciton scattering2L2£ and 
Bose condensation of excitonsi^ 

To compute (\In\ 2 )i3 we define a density matrix that 
contains both radiating and diagonal constraints 

p(RN, Rjv) = Prad(RjV, Rjv) + Pdiag(RjV, Rjv) J (H) 

where 

Pr&d(R-N, R'at) ^(RjV.R'jvl^R-N-l — RjV-l) 



(12) 



and 



Pdiag(RjV, Rjv) — Pn(R<N, Rjv) S(Rn — Rjy) (13) 

Since the radiating and diagonal constraints form two 
disjoint subsets in configuration space we can write the 
probability of being in cither state as 



P(radiating/diagonal state) = 

J % Vrad/diag^iV, Rjv) dRjvdR; 



(14) 



where Z = J /5(Rjv,R^y) dHxdIl' N . Combining Eqs. Q 
and l|14|) . we get an expression suitable for evaluation 
within PIMC: 



{\In\ 2 ) 



P(radiating state) 
P(diagonal state) 



(15) 



In our simulations we use a path integral expansion of 
the density matrix px with a finite number of imaginary 
time slices. The configuration space of this expansion is 



Rat, R^ , . . . , R^" 1 = Rjv) where m is the num- 



ber of time slices. We sample the probability distribution 
Z~ 1 p using the Metropolis algorithm. Since the number 
of time slices m is of order 10 4 in a typical calculation, 
it is essential to use a multilevel Metropolis algorithm^ 
especially when changing the configuration from radiat- 
ing to diagonal state and vice versa. The probability of 
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being in either state can then be estimated from the rel- 
ative frequencies x ra( j and x^iag = 1 — £ r ad of radiating 
and diagonal path configurations in the Markov chain. 

Finally, we arrive at (\In\ 2 )(3 ~ iCrad/^diagi an d from 
Eq. JfJJ) we get the radiative recombination rate, 



(a) 
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2nE gap Epe 
3ft 2 c 3 m 
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When calculating rates, we use the exciton energies from 
the simulation for E gap . 

Since the temperature k^T in our simulations is small 
compared to the single particle level spacing in the dot, 
we can assume that electrons and holes in the biexciton 
are in a singlet state. Therefore, the fcrmion sign problem 
does not occur in our calculations. (For a review on the 
origin of the sign problem, see e.g. Ref. loOh 



B. Configuration Interaction Calculations 

To demonstrate our method, we have also performed 
CI calculations on the same EMA models, Eq. (|TJ. The 
single particle states are calculated by finite-difference 
discretization in a cylindrical cell with 30 nm height 
and 100 nm diameter, with grid spacing of 0.5 nm and 
0.8 in the vertical and radial directions, respectively; 
Coulomb integrals are evaluated by successive over re- 
laxation. This is the same approach used to calculate 
multi-exciton states reported in Ref. l3ll 

For the simulation of excitons and biexcitons in self- 
assembled QDs, our CI expansion uses a 6s5p4d3f2g2hli 
basis set including 44 single particle states. In contrast, 
the CI expansion by Corni et al. uses a 4s4p3d basis set 
(18 single particle states) in the xy-plane and only a sin- 
gle state for the z-direction consisting of Gaussians cen- 
tered on the dot. Contrary to the direct expansion of the 
many-particle wave-function in our approach, Corni et 
al. use this basis set first to solve the restricted Hartree- 
Fock equation to obtain an optimized basis set for the 
CI expansion. Like our approach, the CI calculations 
by Narvaez et al. also use single particle states from a 
non-interacting Hamiltonian to expand the many- body 
wave function. No basis set size is given in Ref. l2ll but 
a previous paper by the same authors using an identical 
method used 6 electron and 10 hole states (12 electron 
and 20 hole states including spin)i22 



IV. TESTS ON PARABOLIC DOT 

To compare our methods, we first consider a model 
system consisting of two oppositely charged particles in 
a harmonic oscillator potential ("Hooke exciton"). The 
Hamiltonian then reads 



H 
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FIG. 2: (Color online) Comparison of results for Hooke ex- 
citon from PIMC and CI (expansion in 54 basis states): (a) 
total and binding energies and (b) |/jv| 2 from integration of 
ODE, from PIMC and CI (dashed line, circles and x's). En- 
ergies are given in units of reduced Hartrees Ha* = /ie 4 /h 2 e 2 , 
where fi is the reduced mass of the exciton. Error bars are of 
the order of symbol size. 



Using center-of-mass and relative coordinates, the prob- 
lem reduces to a ordinary differential equation (ODE) 
that can be integrated numerically with almost arbitrary 
exactness. We apply both the PIMC and CI techniques 
to this system. Within the path integral calculations we 
use a temperature of (3 = 10 Ha* -1 , which is low enough 
to ensure that only the ground state contributes, and 
to = 500 time slices. The CI calculations use 54 single 
particle states to expand the two-particle wave function. 

Figure Ufa) shows the total and the binding energy of 
the two particles. Both CI and PIMC show very good 
agreement with the results from numerical integration of 
the ODE. However, for |ijv| 2 only PIMC shows good con- 
vergence whereas the CI results are in general too low, 
up to a factor of 2 in the weak confinement case. But 
even for strong confinement there is a considerable dis- 
crepancy, although the confinement energy significantly 
exceeds the exciton binding energy (see Fig. H£b)). In 
our calculations we have also found that the CI result 
for | In 1 2 approaches the correct value rather slowly with 
increasing basis set size, thus leading to a false impres- 
sion of convergence. If only the dependence of the result 
on the basis set is used as a measure of convergence, it 
is hard to decide whether a calculation has converged or 
not. 

The true many particle wave function for Coulombic 
interactions must have a coalescence cusp for ri = r 2 r^ 
but a CI expansion of the wave function in products 
of smooth single particle basis functions cannot have a 
non-analytic behaviour. Convergence problems of CI 
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associated with the failure to reproduce this cusp are 
for exam ple solved by using correlated basis functions 
(e.g. Ref. I34l35l36h . It is therefore not surprising that 
CI calculations give better results for energies than for 
the overlap matrix element: The energy is calculated us- 
ing the wave function at every grid point, whereas for 
the decay rate mainly the cusp at ri = r2 enters. That 
the overlap matrix clement is more sensitive to errors 
in the wave function than the energy, also shows in the 
fact that even a Hartree-Fock calculation gets up to 95% 
of the exciton binding energy^ but completely lacks the 
correlation cusp, leading to a decay rate that does not 
depend on the dot size^i Still, it is striking that the CI 
results are able to reproduce energies very accurately and 
yet completely fail to obtain the correct overlap matrix 
element. PIMC does not suffer from a finite basis set and 
can thus reproduce both energies and the overlap matrix 
element very accurately. 

V. RESULTS FOR SELF-ASSEMBLED DOTS 

We have applied these techniques to common single- 
band effective mass models of quantum dots, summarized 
in Table [I] We chose these materials and sizes because 
of availability of published experimental values. The dot 
geometry is a lens shape, with a height to diameter ratio 
of 1:10. The calculations include a wetting layer, modeled 
as a quantum well with thickness i\yL extending from the 
base of the dot. The dot potential consists of potential 
steps of finite height V e (Vh) for electrons (holes) at the 
boundaries of the lens and the wetting layer. The three 
systems we have studied arc: 

1. InGaAs/GaAs: This is the most studied material 
for optical properties of self-assembled dots, and we 
are comparing our results with four separate PL 
rate experiments. Some of these dots are grown 
as alloyed InGaAs material, while others are nom- 
inally pure In As. Even for nominally pure dots, 
intermixing and annealing at high temperatures 
often leads to dots with significant Ga content. 
Based on reported growth conditions and PL en- 
ergies, we have chosen to simulate dots composed 
of Ino.sGao.s As. The dot diameters, from 10 nm to 
60 nm, cover the size range for nearly all dots of this 
material reported in PL studies. We have included 
a 6 monolayer (ML), or 16 A, Ino.5Gao.5As wetting 
layer under the dot pi To show the influence of the 
wetting layer we also give results for i-yvL = 0. 

2. InP/InGaP: We have included a 2 ML, or 5 A, InP 
wetting layer under the dot. 

3. CdSe/ZnSe: We have included a 2 ML, or 5 A, 
CdSc wetting layer under the dot4i 

While our path integral formalism allows for a ther- 
mal distribution of initial states, we have chosen a low 
temperature (T ss 8K) so that we consider only emission 



from the ground state. We have discretized imaginary 
time in the path integral in steps of t = 1.3 x 10~ 5 -£T -1 . 
The simulation time for one dot diameter was approx. 
200 min for the exciton and 350 min for the biexciton on 
10 Athlon MP 1600+ processors. 

In Fig.|21we present results of our path integral calcula- 
tions, along with our CI results and published experimen- 
tal data points. In Fig.^b) we see that the absolute exci- 
ton decay rate Tx increases for large dots with increasing 
dot diameter due to the larger exciton coherence volume. 
As already expected from our model calculations, the CI 
results for the decay rate suffer from underconvergence, 
although the exciton energies from PIMC and CI agree 
very well. This is particularly evident in the decay rate 
for CdSe/ZnSe where the CI result begins to saturate for 
large dot sizes due to missing correlation, whereas the 
Monte Carlo result still increases uniformly with increas- 
ing dot diameter. A similar flattening of the decay ratio 
with increasing dot diameter can also be observed in the 
CI results of Corni et ai, possibly indicating missing cor- 
relation at larger dot sizes. For small dots, we observe 
a minimum of the decay rate when the dot height h be- 
comes comparable to the wetting layer thickness tyvL- 
The InGaAs/GaAs dots without wetting layer do not 
show this behaviour. Note that in this case we only give 
results down to dot diameter d = 15 nm, because the 
exciton becomes unbound for smaller dot sizes. The ex- 
citon decay rate in the InGaAs/GaAs material system 
is larger for tyvL = 16 A than for i-yvL = A since the 
effective dot size and thus the exciton coherence volume 
is larger for the dots including a wetting layer. 

The relative decay rate T xx x of the biexciton, 
Fig. |ija), varies from approx. 2 down to 1.5 for In- 
GaAs/GaAs and InP/InGaP and even down to 1 for the 
CdSe/ZnSe material system. For large dots we observe 
a decrease of Txx/^x for increasing dot size, at small 
dot sizes there is a maximum of the relative biexciton 
decay rate corresponding to the minimum in the exciton 
decay rate. Again, the data for InGaAs/GaAs without 
wetting layer does not show an extremum for small dots 
but reaches Txx/^x ~ 2, corresponding to the strong 
confinement or non- interacting limit. As we explained 
above, CI underestimates decay rates. However, since it 
underestimates both exciton and biexciton decay rates, 
the relative ratio from CI is actually rather similar to the 
Monte Carlo result. 

To gain more insight into the size-dependence of the 
decay rates, it is useful to study the spatial extent of 
the exciton wave function in the dot. The decay rate is 
closely linked to the coherence volume and thus to the 
volume that is filled by the exciton wave function. In 
Fig. 0] we present the size-dependence of the exciton ra- 
dius ax and the standard deviation of the exciton center 
of mass (com) coordinate Ar com = sj (r"* om ) - (r com ) 2 
for the InGaAs/GaAs and CdSe/ZnSe material systems. 
The results for InP/InGaP quantum dots are similar to 
those for InGaAs/GaAs, just as in Fig. EI 
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TABLE I: Parameters used in the calculations. 



dot /barrier 


E h e ^(eV) 


c 


m e 


-11 




AT4 (eV) 


AVfc (eV) 


i\¥L (A) 


E P (eV) 


Ino.5Gao.5As/GaAs 


1.519° 


12. 5 b 


0.067 6 


0.11 6 


0.38 b 


0.250 c 


0.200 c 


16, 


25.7 6 


InP/Ino.sGao.sP 


1.920 d 


12. 6 C 


0.079 c 


0.150 c 


0.600 e 


0A20 d 


0.070 d 


5 


20.4 e 


CdSe/ZnSe 


2.820" 


9.3 ; 


0.130 7 


0.380 9 


1.000 9 


0.735 h 


0.135 h 


5 


17.5 9 



a Ref.|33 

'We approximate the strained InGaAs material in the dot by just 
taking the bulk GaAs value™ 
c Estimated from strain-modified band offsets plotted in Ref . l3St 
d Estimatcd from EPM/VFF calculations (Nair, unpublished). 
e Bulk InP valued 
^BulkCdSe values^ 

h CdSe/ZnSe band offsets chosen to match simulations in Ref. |4Q| . 



For small dot sizes wc find a minimum of both ax and 
Ar com corresponding to the minimum in the decay rate. 
As the dot height becomes comparable to the wetting 
layer thickness, the exciton center of mass moves into 
the wetting layer and the wave function extends further 
into the quantum well underneath the dot. A dot height 
less than £ WL thus corresponds to an effectively larger 
dot size. The increased coherence volume leads to an 
increase in the decay ratio Tx and a decrease of the rel- 
ative biexciton ratio Txx/^x- In the limit of zero dot 
size, we would be in the quantum well situation, how- 
ever, the dipole approximation leading to Eq. JSJ does 
not hold for an extended quantum well stated In the 
case of the QD without wetting layer the coherence vol- 
ume decreases monotonically until the exciton becomes 
unbound, thus we do not observe an extremum in the 
decay rates. 

In the case of larger dot sizes we observe different be- 
haviour for InGaAs/GaAs and CdSe/ZnSe: In the case 
of InGaAs/GaAs, Ar com < ax for all studied diameters 
and the exciton radius ax does not show saturation to- 
wards the bulk value. With respect to these properties of 
the wave function, the dots are in the strong confinement 
limit, whereas the decay rate shows signatures of strong 
to intermediate confinement: The relative biexciton ra- 
tio can be tuned over a relatively large range — from 2 to 
1.5 — by changing the dot geometry, an effect entirely due 
to electronic correlation. 

For CdSe/ZnSe we find a crossover to Ar com > ax 
with increasing dot diameter. This corresponds to the 
weak confinement limit — an exciton "bouncing" around 
in the dot — and thus we observe relative biexciton ratios 
down to 1. However, reported photoluminescence mea- 
surements on this type of QDs arc usually carried out on 
dots with a diameter d around 10 nm, so that the rel- 
evant experimental data for CdSe/ZnSe also lies in the 
strong to intermediate confinement regime. 

When comparing to experiment, we notice that for 
some of the reported data our calculated exciton ener- 
gies are much smaller than the experimental values. In 
these cases (Refs.0@) the growth conditions enhanced al- 
loying. Our parameters for the band offsets do not seem 



to describe these shallow dots very well. However, since 
we use an abrupt potential step for the QD boundaries, 
the step height should not influence our calculated decay 
rates significantly, as long as the exciton is still bound. 
The step height just determines the exponential decay 
length of the wave function into the barrier, so its influ- 
ence on the wave function inside the dot is only indirect. 
Comparison of our calculated results with these experi- 
ments is thus still valid. 

In comparing to the work of Thompson et a/.^fl we 
found it was necessary to re-identify their reported exci- 
ton spectra line as a charged exciton. Our concerns were 
their reported negative (blue-shifted) biexciton binding 
energy and their ratio Txx/^x ~ 2.3. Their spectra 
are very similar to spectra reported by Lomascolo et al., 
which arc dominated by charged exciton^ only with the 
exciton/charged exciton labels switched. Since Thomp- 
son et al. called their identification of the charged and 
neutral exciton tentative, and did not offer any alterna- 
tive explanation for the unusual energy shift and relative 
decay rate of their supposed neutral exciton/biexciton 
pair, we chose to compare to the data they had attributed 
to the charged exciton and charged biexciton. These 
states have Txx/^x ~ 2.06 and a biexciton binding en- 
ergy of +2 mcV. This identification makes little difference 
in our comparison of decay rates, but does give us much 
better agreement for the positive biexciton binding en- 
ergy and is consistent with the relative biexciton decay 
rate in the strong-confinement limit. Thompson et al. do 
not give any value for the dot size in their experiment, but 
claim that their dots are smaller than those of other pho- 
toluminescence experiments. We chose to attribute their 
data a dot diameter of 20 nm, based on our energy calcu- 
lations for the InGaAs/GaAs dot without wetting layer. 
In the absence of further information about the wetting 
layer thickness in the experiment, this seems reasonable 
since the nominal InAs coverage in this experiment is 
only 1.7ML. Accordingly, we will always compare the ex- 
perimental data from Thompson et al. as well as the data 
from Ulrich et ali^ (experimental wetting layer thickness 
1 ML) with the results of our calculation with i-yvL = 
ML. 
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FIG. 3: (Color online) Our results of path-integral calculations for for InGaAs/GaAs, InP/InGaP, and CdSe/ZnSe lens-shaped 
dots of different diameters, with height/diameter = 0.1. Rows of panels present (a) the relative decay rate of biexciton to 
exciton, (b) the absolute decay rate of exciton, (c) the exciton energy, at which the exciton luminescence peak would be 
observed, and (d) the energy shift of the biexciton luminescence peak. Solid circles are path integral results for twh > 0, open 
circles for 4wl = 0. Error bars are only given if the error exceeds the symbol size. We compare to our configuration interaction 
calculations for 4wl > (x's), which miss some of the correlation. We also show experimental data for InGaAs/GaAs (A's 
are data from Ref. llOl where we have assumed d = 20 and used their measured X* and XX* data points; CPs are data from 
Ref. U 0's are from Ref. 111! and +'s are from Ref. Il2t) and for CdSe (CPs are data from Ref. 0, the range of data from Ref. 
is indicated by dashed crosses) . 



Our PIMC results for the exciton decay rate in In- 
GaAs/GaAs agree with experiment within about a fac- 
tor of two, but seem to systematically overestimate the 
decay rate. This could be due to our simplified model 
of an ideal dot. The agreement with experiment could, 
for example, be improved by including the effects of al- 
loying in the dot potential. Disorder introduced by al- 
loying leads to stronger localization of the particles in 
the dot^ thus reducing the coherence volume and the 
electron- hole overlap. Another possiblity for improving 
the path integral results would be to use a dot potential 
from strain calculations^ since strain might also lead to 



an increased electron-hole separation^* The inclusion of 
such single-particle potentials in PIMC is perfectly feasi- 
ble and does not introduce any additional computational 
cost. Even in the strong confinement or non-interacting 
limit | In | 2 ~ 1, so Tx ~ 2 ns _1 using the parameters 
from Table |U Thus the low decay rate from Refs. ITol and 
[Til cannot be explained in our model, hinting at the need 
of a more detailed dot potential. However, for the study 
of the size dependence of the decay rates a model poten- 
tial is perfectly valid and yields results that are easier to 
interpret. 
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FIG. 4: (Color online) Exciton radius ax (D's) and exciton 
center of mass fluctuations Ar com (A's) for InGaAs/GaAs 
dots with a 6 ML wetting layer (full symbols) and for 
CdSe/ZnSe dots (open symbols). 



The Monte Carlo calculations can reproduce the range 
of observed relative biexciton ratios from 2 to 1.5. The 
data from Refs. ^ and ITU is described very well by the 
QD without wetting layer, whereas the data from Ref. 
seems to be best reproduced by a QD with wetting layer, 
although we have to assume a somewhat larger effec- 
tive dot size. The calculated biexciton binding ener- 
gies are also close to the experimental values. The ex- 
tremely low relative biexciton decay rate from Ref. IT2I 
T xx /T x = 0.33 however cannot be explained at all in 
our model. In the original paper, the low biexciton decay 
rate was attributed to weak confinement effects, but from 
our calculations we can conclude that InGaAs/GaAs QDs 
with diameters around 50 nm are still far from the weak 
confinement regime. 

We are not aware of any studies on the exciton and 
biexciton dynamics in single InP/InGaP QDs, but our 
calculations arc within the reported exciton life time 
range of 100-500 ps for QD ensembles with dot diam- 
eters between 20 nm and 40 nmi^ 

When comparing to the experimental data on single 
CdSe/ZnSe dots by Patton et al.jL we chose to only give 
the range of the reported data. The variation of exci- 
ton energies in Patton et al. is attributed to different 
localization potentials, i.e. different Cd concentration, 
and not to different dot sizes. The dot diameter for 
their samples is reported to be between 5 nm and 10 
nm4£ The CdSe/ZnSe exciton decay rates calculated by 
PIMC agree very well with the reported experimental 
data. However, we completely fail to reproduce the rel- 
ative biexciton ratio Txx/Tx ~ 1 from Bacher et al& 
Such a low biexciton decay rate is only to be expected 
for very large dots in our simulation. The experiments 
by Patton et al. on very similar sized QDs in contrast 
yielded a relative biexciton ratio Txx /T x ~ 2 with a 



rather large experimental spread (Txx /Fx ~ 1.4 — 2.8). 
We expect CdSe/ZnSe QDs of about 10 nm to be towards 
the strong confinement limit, consistent with the exper- 
iment by Patton et al. Therefore we presume that more 
knowledge about the dot potential would be needed to 
explain the results of Bacher et al., a simple box model, 
as suggested in Ref. 0, is certainly not enough. The fact 
that the exciton energies from Ref. are well explained by 
our model whereas the rather shallow dots from Ref. are 
not, supports this presumption. It should also be noted 
that the QDs of Ref. lisl were grown under conditions sim- 
ilar to those of Bacher et al. They show recombination 
energies and exciton lifetimes comparable to the results 
of Bacher et at, but much shorter biexciton lifetimes. 
Since no estimate of the dot size was given, wc cannot 
directly compare to our calculations, but the reported ra- 
tio of Txx /Tx ~ 1.4 agrees well with the ratios expected 
from our calculation. 

We cannot compare our hitherto obtained rates with 
the calculations of Corni et al. on the sizc-dcpcndcncc 
of the exciton and biexciton decay rates because of the 
different dot potentials. However, if wc apply our tech- 
nique to the truncated parabola potentials used in their 
study, we obtain exciton decay rates that are for small 
dots around 50%, and for large dots even up to two times 
larger than the results of Ref. |2(J Given that even our 
CI expansion, using a large basis set of 44 single particle 
states, yields absolute rates that are too low, it is not sur- 
prising that the much smaller basis set of Corni et al. also 
fails to calculate absolute rates. Yet, the CI calculations 
show the right trends — increasing decay rate and decreas- 
ing Txx/Tx with increasing dot size — compatible with 
our results. Also, the relative biexciton decay rate from 
PIMC is very similar to the one obtained by Corni et al. 

Narvaez et al. have performed CI calculations on the 
height dependence of recombination rates in lens-shaped 
Ino.6Gao.4As/GaAs quantum dots with a fixed diameter 
(25.2 nm), using atomic pseudopotentials and a realistic 
model for alloying. Their calculates decay rates lie in the 
range of 0.4-0.5 ns _1 , a factor of four lower than our re- 
sults, but also a factor of two lower than experimental 
valucs^SiAi Their reported relative biexciton decay ratio 
Txx/Tx — 4 is a factor of 2 larger than what is ex- 
pected for strong confinement and is to our knowledge 
not observed in experiment. Path integral techniques 
cannot be adapted to using pseudopotentials easily, and 
thus we cannot directly compare results. Still, from our 
experience with the calculation of rates from CI, we are 
somewhat concerned with the absolute value of the rates. 
The reported basis set size in Ref. |3^ is much smaller than 
the one used in this study. For example, their minimum 
of the exciton lifetime (corresponding to a maximum in 
the decay rate) at a dot height of 65A could possibly 
be a sign of missed correlation at larger dot sizes: From 
our calculations we would expect the decay rate to grow 
monotonically with dot volume. The decrease of the ex- 
citon lifetime found at smaller dot heights however is 
compatible with our findings. 
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VI. CONCLUSION 

We have developed a path integral Monte Carlo ap- 
proach for studying exciton and biexciton recombination 
rates in self-assembled quantum dots. This technique 
allows us to study general 3-d potentials for a wide va- 
riety of single-band EMA models. Our calculations in- 
dicate that self-assembled dots are in the strong to in- 
termediate confinement regime, where Coulomb correla- 
tion effects are becoming important. In particular, for 
large dots we see a clear monotonic rise in recombina- 
tion rate versus diameter, and a decrease in the relative 
biexciton decay rate, Txx/^x- From our calculations 
we can state that relative decay rates Txx /Tx ~ 1.5 — 2 
are expected for typical photoluminescence experiments, 
an effect due entirely to correlation. We have seen that 
quantum dots of the size used in PL experiments tend 
towards the strong confinement regime. Thus, the low 
relative biexciton decay rates Txx /Tx < 1 from some 
experimentsjSii^ that were attributed to weak confine- 
ment, cannot be explained by weak confinement effects. 
It should be noted that in single dot experiments rather 
large dot-to-dot fluctuations have been reported^ Given 
the spread of experimental data, our calculations com- 
pare rather favorably against experiment. 

We have further shown that CI expansions using un- 
corrected single particle basis sets have severe shortcom- 
ings in calculating decay rates. Rather surprisingly, we 
have found that CI expansion in a large basis set of 44 
states underestimate decay rates by far, even for dot sizes 



comparable to the exciton Bohr radius and although the 
calculated energies were well-converged. Yet, due to a 
cancellation of errors, the relative biexciton decay rate 
calculated by CI was found to be similar to the path 
integral result. Also, trends were in general reproduced 
correctly by CI. CI has some advantages over PIMC, such 
as being able to use an atomic description of the quan- 
tum dot. However, absolute decay rates from CI must be 
regarded with caution. 

In conclusion, we have developed a microscopic path- 
integral technique for calculating exciton and biexciton 
decay rates, that fully treats quantum correlation in re- 
alistic models. We can apply this to arbitrary geome- 
tries within single-band EMA models. Our calculations 
on lens-shaped self-assembled dots indicate that these 
commonly studied structures are in the regime between 
strong and intermediate confinement. The formalism has 
a built-in thermal distribution of carriers that we have 
not yet exploited. Another area for future research is an 
extension of this technique to semiconductors with indi- 
rect band gaps, such as Si/Ge. 
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